library(tidyverse)
library(rlang)
library(leaflet)
hospitals <- read_csv("../data/hospital_list.csv")
Parsed with column specification:
cols(
.default = col_character(),
zip = [32mcol_double()[39m,
`Phone Number` = [32mcol_double()[39m,
score = [32mcol_double()[39m,
longitude = [32mcol_double()[39m,
latitude = [32mcol_double()[39m,
population = [32mcol_double()[39m,
poverty = [32mcol_double()[39m,
percent = [32mcol_double()[39m,
poverty_percent = [32mcol_double()[39m
)
See spec(...) for full column specifications.
hospitals
hospital_locations <- hospitals %>%
select(state, longitude, latitude)
write_rds(hospital_locations, "../data/hospital_locations.rds")
hospital_locations
population <- usmap::countypop
hospital_count <- hospitals %>%
count(state, original_county)
state_hospitals <- population %>%
left_join(hospital_count, by = c("abbr" = "state", "county" = "original_county")) %>%
mutate(hospitals = ifelse(is.na(n), 0, n)) %>%
select(
fips,
state = abbr,
county,
population = pop_2015,
hospitals
)
state_hospitals
set.seed(100)
hospital_sample <- state_hospitals %>%
sample_frac(0.3)
hospital_sample
model <- lm(hospitals ~ population, data = hospital_sample)
write_rds(model, "../data/model.rds")
summary(model)
Call:
lm(formula = hospitals ~ population, data = hospital_sample)
Residuals:
Min 1Q Median 3Q Max
-10.7909 -0.7106 0.0908 0.2498 8.3805
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.027e-01 3.942e-02 17.83 <2e-16 ***
population 7.907e-06 1.182e-07 66.90 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.149 on 941 degrees of freedom
Multiple R-squared: 0.8263, Adjusted R-squared: 0.8261
F-statistic: 4475 on 1 and 941 DF, p-value: < 2.2e-16
predictions <- predict(model,
newdata = state_hospitals,
interval = "prediction") %>%
as_tibble() %>%
mutate_all(round)
predictions
hospital_results <- state_hospitals %>%
bind_cols(predictions) %>%
mutate(result = case_when(
hospitals < lwr ~ -1,
hospitals > upr ~ 1,
TRUE ~ 0)) %>%
mutate(
county_key = str_remove_all(county, " County"),
county_key = str_remove_all(county_key, " Parish"),
county_key = str_remove_all(county_key, " city"),
county_key = str_to_lower(county_key),
county_key = str_remove_all(county_key, " "),
county_key = str_replace_all(county_key, "st. ", "saint"),
)
write_rds(hospital_results, "../data/hospitals.rds")
hospital_results
State map
shapes <- read_rds("../data/shapes.rds") %>%
mutate(
county_key = str_to_lower(county),
county_key = str_remove_all(county_key, " "),
county_key = str_replace_all(county_key, "st. ", "saint")
)
shapes
county_hospitals <- hospital_results %>%
select(- county) %>%
left_join(shapes, by = c("state", "county_key")) %>%
filter(state != "PR", state != "AK") %>%
filter(!is.na(county))
county_hospitals %>%
mutate(popup = paste0("County:", county, " Population:", population, " Hospitals:", hospitals)) %>%
select(popup, everything())
states <- county_hospitals %>%
group_by(state) %>%
summarise() %>%
pull()
under <- "#CC79A7"
over <- "#0072B2"
at_level <- "#008b00"
hospital_color <- "#F0E442"
st <- c("LA")
#st <- states
include <- c(-1, 1, 0)
add_hospitals <- 0
initial_map <- leaflet() %>%
addProviderTiles(providers$CartoDB) %>%
addLegend("bottomright",
color = c(under,over ,at_level, hospital_color),
labels = c("Less hopitals than expected",
"More hospitals than expected", "Within Range",
"Hospital Location"),
title = "Legend",opacity = 0.5)
locations <- hospital_locations %>%
filter(state %in% st) %>%
select(longitude, latitude) %>%
mutate_all(round, 3) %>%
count(longitude, latitude) %>%
mutate(popup = paste0("Hospitals: ", n))
counties <- county_hospitals %>%
filter(state %in% st,
result %in% include) %>%
mutate(popup = paste0("<b>", county, "</b>",
"<br/>Population: ", prettyNum(population, big.mark = ","),
"<br/>Hospitals: ", hospitals,
"<br/>Expected: ", fit
)) %>%
mutate(color = case_when(
result == 0 ~ at_level,
result == 1 ~ over,
result == -1 ~ under
))
county_map <- counties %>%
transpose() %>%
map(~{
county <- .x
transpose(county$data) %>%
map(~list(
data = .x$data,
color = county$color,
popup = county$popup
))
}) %>%
flatten() %>%
reduce(
~ addPolygons(.x,
lng = .y$data$long,
lat = .y$data$lat,
color = .y$color,
popup = .y$popup,
weight = 1, fillOpacity = 0.3),
.init = initial_map)
if(add_hospitals == 1) {
county_map <- county_map %>%
addCircleMarkers(
lng = locations$longitude,
lat = locations$latitude,
radius = 1.5 *locations$n,
popup = locations$popup,
fillColor = hospital_color, color = "gray",
fillOpacity = 0.7,weight = 1)
}
county_map
NA
LS0tDQp0aXRsZTogIkFjY2VzcyB0byBjYXJlIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KYGBge3J9DQpsaWJyYXJ5KHRpZHl2ZXJzZSkNCmxpYnJhcnkocmxhbmcpDQpsaWJyYXJ5KGxlYWZsZXQpDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbHMgPC0gcmVhZF9jc3YoIi4uL2RhdGEvaG9zcGl0YWxfbGlzdC5jc3YiKQ0KaG9zcGl0YWxzDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbF9sb2NhdGlvbnMgPC0gaG9zcGl0YWxzICAlPiUNCiAgc2VsZWN0KHN0YXRlLCBsb25naXR1ZGUsIGxhdGl0dWRlKQ0Kd3JpdGVfcmRzKGhvc3BpdGFsX2xvY2F0aW9ucywgIi4uL2RhdGEvaG9zcGl0YWxfbG9jYXRpb25zLnJkcyIpDQpob3NwaXRhbF9sb2NhdGlvbnMNCmBgYA0KDQpgYGB7cn0NCnBvcHVsYXRpb24gPC0gdXNtYXA6OmNvdW50eXBvcA0KDQpob3NwaXRhbF9jb3VudCA8LSBob3NwaXRhbHMgJT4lDQogIGNvdW50KHN0YXRlLCBvcmlnaW5hbF9jb3VudHkpIA0KDQoNCnN0YXRlX2hvc3BpdGFscyA8LSBwb3B1bGF0aW9uICU+JQ0KICBsZWZ0X2pvaW4oaG9zcGl0YWxfY291bnQsIGJ5ID0gYygiYWJiciIgPSAic3RhdGUiLCAiY291bnR5IiA9ICJvcmlnaW5hbF9jb3VudHkiKSkgJT4lDQogIG11dGF0ZShob3NwaXRhbHMgPSBpZmVsc2UoaXMubmEobiksIDAsIG4pKSAlPiUNCiAgc2VsZWN0KA0KICAgIGZpcHMsIA0KICAgIHN0YXRlID0gYWJiciwgDQogICAgY291bnR5LCANCiAgICBwb3B1bGF0aW9uID0gcG9wXzIwMTUsDQogICAgaG9zcGl0YWxzDQogICAgKQ0KDQpzdGF0ZV9ob3NwaXRhbHMNCmBgYA0KDQoNCg0KDQpgYGB7cn0NCnNldC5zZWVkKDEwMCkNCg0KaG9zcGl0YWxfc2FtcGxlIDwtIHN0YXRlX2hvc3BpdGFscyAlPiUNCiAgc2FtcGxlX2ZyYWMoMC4zKQ0KDQpob3NwaXRhbF9zYW1wbGUNCmBgYA0KDQpgYGB7cn0NCm1vZGVsIDwtIGxtKGhvc3BpdGFscyB+IHBvcHVsYXRpb24sIGRhdGEgPSBob3NwaXRhbF9zYW1wbGUpDQp3cml0ZV9yZHMobW9kZWwsICIuLi9kYXRhL21vZGVsLnJkcyIpDQpzdW1tYXJ5KG1vZGVsKQ0KYGBgDQoNCmBgYHtyfQ0KcHJlZGljdGlvbnMgPC0gcHJlZGljdChtb2RlbCwgDQogIG5ld2RhdGEgPSBzdGF0ZV9ob3NwaXRhbHMsIA0KICBpbnRlcnZhbCA9ICJwcmVkaWN0aW9uIikgJT4lDQogIGFzX3RpYmJsZSgpICU+JQ0KICBtdXRhdGVfYWxsKHJvdW5kKQ0KDQpwcmVkaWN0aW9ucw0KYGBgDQoNCmBgYHtyfQ0KaG9zcGl0YWxfcmVzdWx0cyA8LSBzdGF0ZV9ob3NwaXRhbHMgJT4lDQogIGJpbmRfY29scyhwcmVkaWN0aW9ucykgJT4lDQogIG11dGF0ZShyZXN1bHQgPSBjYXNlX3doZW4oDQogICAgaG9zcGl0YWxzIDwgbHdyIH4gLTEsDQogICAgaG9zcGl0YWxzID4gdXByIH4gMSwNCiAgICBUUlVFIH4gMCkpICU+JQ0KICBtdXRhdGUoDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eSwgIiBDb3VudHkiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBQYXJpc2giKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBjaXR5IiksDQogICAgY291bnR5X2tleSA9IHN0cl90b19sb3dlcihjb3VudHlfa2V5KSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiAiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlcGxhY2VfYWxsKGNvdW50eV9rZXksICJzdC4gIiwgInNhaW50IiksDQogICAgICApICANCg0Kd3JpdGVfcmRzKGhvc3BpdGFsX3Jlc3VsdHMsICIuLi9kYXRhL2hvc3BpdGFscy5yZHMiKQ0KDQpob3NwaXRhbF9yZXN1bHRzIA0KYGBgDQoNCiMjIFN0YXRlIG1hcA0KDQpgYGB7cn0NCnNoYXBlcyA8LSByZWFkX3JkcygiLi4vZGF0YS9zaGFwZXMucmRzIikgJT4lDQogIG11dGF0ZSgNCiAgICBjb3VudHlfa2V5ID0gc3RyX3RvX2xvd2VyKGNvdW50eSksDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eV9rZXksICIgIiksDQogICAgY291bnR5X2tleSA9IHN0cl9yZXBsYWNlX2FsbChjb3VudHlfa2V5LCAic3QuICIsICJzYWludCIpDQogICkgDQoNCnNoYXBlcyANCmBgYA0KYGBge3J9DQpjb3VudHlfaG9zcGl0YWxzIDwtIGhvc3BpdGFsX3Jlc3VsdHMgICU+JQ0KICBzZWxlY3QoLSBjb3VudHkpICU+JQ0KICBsZWZ0X2pvaW4oc2hhcGVzLCBieSA9IGMoInN0YXRlIiwgImNvdW50eV9rZXkiKSkgJT4lDQogIGZpbHRlcihzdGF0ZSAhPSAiUFIiLCBzdGF0ZSAhPSAiQUsiKSAlPiUNCiAgZmlsdGVyKCFpcy5uYShjb3VudHkpKQ0KDQpjb3VudHlfaG9zcGl0YWxzICU+JQ0KICBtdXRhdGUocG9wdXAgPSBwYXN0ZTAoIkNvdW50eToiLCBjb3VudHksICIgUG9wdWxhdGlvbjoiLCBwb3B1bGF0aW9uLCAiIEhvc3BpdGFsczoiLCBob3NwaXRhbHMpKSAlPiUNCiAgc2VsZWN0KHBvcHVwLCBldmVyeXRoaW5nKCkpDQpgYGANCg0KYGBge3J9DQpzdGF0ZXMgPC0gY291bnR5X2hvc3BpdGFscyAlPiUNCiAgZ3JvdXBfYnkoc3RhdGUpICU+JQ0KICBzdW1tYXJpc2UoKSAlPiUNCiAgcHVsbCgpDQoNCmBgYA0KDQoNCmBgYHtyfQ0KDQp1bmRlciA8LSAiI0NDNzlBNyINCm92ZXIgPC0gIiMwMDcyQjIiDQphdF9sZXZlbCA8LSAiIzAwOGIwMCINCmhvc3BpdGFsX2NvbG9yIDwtICIjRjBFNDQyIg0Kc3QgPC0gYygiTEEiKQ0KI3N0IDwtIHN0YXRlcw0KaW5jbHVkZSA8LSBjKC0xLCAxLCAwKQ0KYWRkX2hvc3BpdGFscyA8LSAwDQoNCmluaXRpYWxfbWFwIDwtIGxlYWZsZXQoKSAlPiUNCiAgYWRkUHJvdmlkZXJUaWxlcyhwcm92aWRlcnMkQ2FydG9EQikgJT4lDQogIGFkZExlZ2VuZCgiYm90dG9tcmlnaHQiLCANCiAgICAgICAgICAgIGNvbG9yICA9IGModW5kZXIsb3ZlciAsYXRfbGV2ZWwsIGhvc3BpdGFsX2NvbG9yKSwgDQogICAgICAgICAgICBsYWJlbHMgPSBjKCJMZXNzIGhvcGl0YWxzIHRoYW4gZXhwZWN0ZWQiLA0KICAgICAgICAgICAgICAgICAgICAgICAiTW9yZSBob3NwaXRhbHMgdGhhbiBleHBlY3RlZCIsICJXaXRoaW4gUmFuZ2UiLCANCiAgICAgICAgICAgICAgICAgICAgICAgIkhvc3BpdGFsIExvY2F0aW9uIiksIA0KICAgICAgICAgICAgdGl0bGUgID0gIkxlZ2VuZCIsb3BhY2l0eSA9IDAuNSkNCg0KbG9jYXRpb25zIDwtIGhvc3BpdGFsX2xvY2F0aW9ucyAlPiUNCiAgZmlsdGVyKHN0YXRlICVpbiUgc3QpICU+JQ0KICBzZWxlY3QobG9uZ2l0dWRlLCBsYXRpdHVkZSkgJT4lDQogIG11dGF0ZV9hbGwocm91bmQsIDMpICU+JQ0KICBjb3VudChsb25naXR1ZGUsIGxhdGl0dWRlKSAlPiUNCiAgbXV0YXRlKHBvcHVwID0gcGFzdGUwKCJIb3NwaXRhbHM6ICIsIG4pKQ0KDQpjb3VudGllcyA8LSBjb3VudHlfaG9zcGl0YWxzICU+JSANCiAgZmlsdGVyKHN0YXRlICVpbiUgc3QsDQogICAgICAgICByZXN1bHQgJWluJSBpbmNsdWRlKSAlPiUNCiAgbXV0YXRlKHBvcHVwID0gcGFzdGUwKCI8Yj4iLCBjb3VudHksICI8L2I+IiwNCiAgICAgICAgICAgICAgICAgICAgICAgICI8YnIvPlBvcHVsYXRpb246ICIsIHByZXR0eU51bShwb3B1bGF0aW9uLCBiaWcubWFyayA9ICIsIiksIA0KICAgICAgICAgICAgICAgICAgICAgICAgIjxici8+SG9zcGl0YWxzOiAiLCBob3NwaXRhbHMsDQogICAgICAgICAgICAgICAgICAgICAgICAiPGJyLz5FeHBlY3RlZDogIiwgZml0DQogICAgICAgICAgICAgICAgICAgICAgICApKSAlPiUgIA0KICBtdXRhdGUoY29sb3IgPSBjYXNlX3doZW4oDQogICAgcmVzdWx0ID09ICAwICB+IGF0X2xldmVsLA0KICAgIHJlc3VsdCA9PSAgMSAgfiBvdmVyLA0KICAgIHJlc3VsdCA9PSAtMSAgfiB1bmRlcg0KICApKSANCg0KY291bnR5X21hcCA8LSBjb3VudGllcyAlPiUNCiAgdHJhbnNwb3NlKCkgJT4lDQogIG1hcCh+ew0KICAgIGNvdW50eSA8LSAueA0KICAgICB0cmFuc3Bvc2UoY291bnR5JGRhdGEpICAlPiUNCiAgICAgICBtYXAofmxpc3QoDQogICAgICAgICBkYXRhID0gLngkZGF0YSwgDQogICAgICAgICBjb2xvciA9IGNvdW50eSRjb2xvciwNCiAgICAgICAgIHBvcHVwID0gY291bnR5JHBvcHVwDQogICAgICAgICApKQ0KICAgIH0pICU+JQ0KICBmbGF0dGVuKCkgJT4lDQogIHJlZHVjZSgNCiAgICB+IGFkZFBvbHlnb25zKC54LCANCiAgICAgICAgICAgICAgICAgIGxuZyA9IC55JGRhdGEkbG9uZywgDQogICAgICAgICAgICAgICAgICBsYXQgPSAueSRkYXRhJGxhdCwgDQogICAgICAgICAgICAgICAgICBjb2xvciA9IC55JGNvbG9yLCANCiAgICAgICAgICAgICAgICAgIHBvcHVwID0gLnkkcG9wdXAsDQogICAgICAgICAgICAgICAgICB3ZWlnaHQgPSAxLCBmaWxsT3BhY2l0eSA9IDAuMyksIA0KICAgIC5pbml0ID0gaW5pdGlhbF9tYXApICANCg0KDQppZihhZGRfaG9zcGl0YWxzID09IDEpIHsNCiAgY291bnR5X21hcCA8LSBjb3VudHlfbWFwICU+JQ0KICAgIGFkZENpcmNsZU1hcmtlcnMoDQogICAgICBsbmcgPSBsb2NhdGlvbnMkbG9uZ2l0dWRlLCANCiAgICAgIGxhdCA9IGxvY2F0aW9ucyRsYXRpdHVkZSwNCiAgICAgIHJhZGl1cyA9IDEuNSAqbG9jYXRpb25zJG4sDQogICAgICBwb3B1cCA9IGxvY2F0aW9ucyRwb3B1cCwNCiAgICAgIGZpbGxDb2xvciA9IGhvc3BpdGFsX2NvbG9yLCBjb2xvciA9ICJncmF5IiwgDQogICAgICBmaWxsT3BhY2l0eSA9IDAuNyx3ZWlnaHQgPSAxKQ0KICB9IA0KDQpjb3VudHlfbWFwDQogIA0KYGBgDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg==